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We examine the thermal fluctuations of the local electric field and the dipole moment //j. 
in liquid water at T = 298 K between metal walls in electric held applied in the perpendicular 
direction. We use analytic theory and molecular dynamics simulation. In this situation, there is 
a global electrostatic coupling between the surface charges on the walls and the polarization in 
the bulk. Then, the correlation function of the polarization density Pz{r) along the applied field 
contains a homogeneous part inversely proportional to the cell volume V. Accounting for the long- 
range dipolar interaction, we derive the Kirkwood-Frohlich formula for the polarization fluctuations 
when the specimen volume v is much smaller than V. However, for not small v/V, the homogeneous 
part comes into play in dielectric relations. We also calculate the distribution of E^ff^ in applied field. 

As a unique feature of water, its magnitude obeys a Gaussian distribution with a large mean 

value Eq = 17 V/nm, which arises mainly from the surrounding hydrogen-bonded molecules. Since 
\pii.\Eo ~ SOksT, fj,). becomes mostly parallel to E^f^. As a result, the orientation distributions 
of these two vectors nearly coincide, assuming the classical exponential form. In dynamics, the 
component of parallel to E]f^{t) changes on the timescale of the hydrogen bonds ~ 5 ps, while 
its smaller perpendicular component undergoes librational motions on timescales of 0.01 ps. 


I. INTRODUCTION 

The dielectric properties of highly polar fluids are com¬ 
plicated because of the long-range dipolar interaction. A 
great number of papers have been devoted on the dielec¬ 
tric constant £ in such fluidJi^. Lorentz and Debye pre¬ 
sented classical formulas for £ using the concept of inter¬ 
nal or local electric field acting on each dipole. However, 
their formulas were inapplicable to highly polar fluids, 
because they did not properly account for the electro¬ 
static interaction between a dipole considered and sur¬ 
rounding ones. Onsager!^ presented a mean-field theory, 
where each dipole is under influence of the reaction field 
produced by a surrounding continuum. KirkwoocP fur¬ 
ther took into account the short-range orientation corre¬ 
lation introducing the so-called correlation factor gx- We 
mention subsequent statistical-mechanical theorie^^^^iil. 

In molecular dynamics simulation^^lt^, e has been cal¬ 
culated on the basis of its linear response expressions un¬ 
der the three-dimensional (3D) periodic boundary condi¬ 
tion. In such studies, the Ewald method has been used 
to efficiently sum the electrostatic interaction 
thermore, we mention simulations on dipoles (and ions) 
in applied electric flelcPHEI] Some group^^MlZI intro¬ 
duced image charges outside the cell to realize constant 
electric potentials at z = 0 and H under the periodic 
boundary condition in the x and y axes. In this paper, 
we use this method. As simplified schemes of applying 
electric field, Yeh and Berkowitd^ assumed empty slabs 
(instead of metal walls) in contact with water neglecting 
the image charges, while Petersen et accounted for 
the primary image charges closest to the walls. On the 
other hand, Willard et used Siepmann and Sprik’s 
electrode modeP^, where atomic particles form a crystal 
in contact with water molecules and their charges vary 


continuously to realize the metallic boundary condition. 

In this paper, we examine the pair correlation func¬ 
tion Gai 3 {ri,r 2 ) = {Pa{'i'i)pp{r 2 )) for the polarization 
density Pair) {a,j3 = x,y,z). For infinite systems, 
FelderhoP found that Gap consists of short-range and 
dipolar parts and its integral within a sphere yields the 
Kirkwood-Frohlich formula. In our theory, there arises a 
global electrostatic coupling between the fluctuations of 
the charge density on the metal walls and those of the 
bulk polarization perpendicular to the walls (along the z 
axis). As a result, the zz component Gzz contains a ho¬ 
mogeneous part independent of ri and r 2 due to the pres¬ 
ence of metal walls. This homogeneous part is inversely 
proportional to the cell volume V and largely alters the 
fluctuations of the total polarization = f drpz{r) 
along the applied field. In our simulation, we also demon¬ 
strate the presence of the angle-dependent dipolar part 
decaying as \ri - r 2 \~^ in Ga 0 {ri,r 2 ). 

The local electric field acting on a point dipole is the 
sum of the contributions from the other dipoles and the 
surface charges on the electrodes. For polar fluids, the 
local field on a molecule k may be defined as follows. 
We change the molecular dipole /Xj, by dpLf. and express 
the resultant change in the electrostatic energy as E^^^ ■ 
dfij.. In our recent paper on watei^, we determined 
in this manner, which is a linear combination of the 
microscopic fields at the three molecular charge sites. In 
this paper, we investigate the microscopic basis of the 
dielectric properties of water with this expression. 

In the classical pictur^, the local field is proportional 
to the applied held and is small. However, in liquid wa¬ 
ter, the microscopic flelds on the molecular charges and 
the molecular local held are all very lar ge being of or¬ 
der 15 V/nm (~ e/cr^ with cr = 3.2 ApUMEl As a 
unique feature of water, such strong electric flelds are 
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mostly produced by the molecules hydrogen-bonded to 
each molecule k. In this paper, we shall see that the 
dipole moment is nearly parallel to the local field 
for each molecule even in applied field. In fact, if we in¬ 
crease the angle between fij. and from 0 to tt/ 2 , 

an energy of order ^ SOfcBT is needed even at 

T = 298 K, where /xq = |/Xj.|. Furthermore, our find¬ 
ing in applied field is that the distribution of the dipole 
angle 0 with respect to the applied field is of the expo¬ 
nential form oc exp(/ioF'p cos ^//cbT') even in the nonlin¬ 
ear regime, where the effective field Ep yields the average 
polarization in accord with the Kirkwood-Frbhlich theory 
in the linear regime. Note that this form was assumed in 
the classical theorjff^t^, but it needs to be justified par¬ 
ticularly for polar fluids with strong local fields. 

In dynamics, we shall see that each dipole /Xj.(t) con¬ 
sists of a slow part on a timescale of 5 ps and a much 
smaller, fast part on a timescale of 0.01 ps. The slow 
part is parallel to and moves with E^^^ {t) , while the fast 
part is perpendicular to E'^'^{t) and undergoes librational 
motions. This indicates that the dipole motions are gov¬ 
erned by the hydrogen bond dynamics i n liquid water. 
In the previous simulations on liquid watepiSIiSESlHi]^ 
tinctly separated slow and fast molecular motions have 
been observed due to the hydrogen bond network. 

The organization of this paper is as follows. In Sec. 
II, we will summarize the theoretical background and 
present a theory on the polarization correlation function 
with electrodes. In Sec. Ill, we will give numerical re¬ 
sults on the thermal fluctuations of the local field and the 
dipole moments, the orientation time-correlation func¬ 
tions, and the polarization correlations. In Appendix A, 
we will give an expression for the local electric field. In 
Appendix B, we will present a continuum theory of the 
polarization correlations to derive the homogeneous part. 
In Appendix C, we will calculate the polarization fluctu¬ 
ations in a spherical specimen under the influence of the 
reaction field. In Appendix D, we will calculate the self¬ 
part of the polarization correlation function for water and 
the radius-dependent Kirkwood correlation factor. 


II. THEORETICAL BACKGROUND 

In this section, we explain our simulation method and 
present analytic theory. Supplementary numerical results 
are also given. 


A. Electrostatics between metal walls 

As in our previous worlP^^ use the TIP4P/2005 
modeP^ without ions, where the number of the water 
molecules is = 2400. The temperature T is fixed at 298 
K in the NVT ensemble with a Nose-Hoover thermostat. 

For each molecule k, its oxygen atom and two pro¬ 
tons are at Vf^o, and respectively. Partial 

charges are at the protons and another point with 


qh = 0.5564e and = —2gH, respectively. The dipole 
moment of molecule k is written as 


Mfe = 9H(^fcHl + T’fcH2 — 2rfeM) = Mo’T'fe, (1) 

where is the unit vector along /x^, with /xq = 2.305D. 
The charge position r^M is slightly shifted from r^o along 
rtk- There is no induced dipole moment in this model. 
See Appendix A for more details. Hereafter, k and £ 
represent molecules, while i = (fc, a) and j = (£, /3) stand 
for molecular charges with a, /3 = HI, H2, M. 

In our system, N molecules are in a Lx Lx H cell with 
L = 41.SA and H = 44.tA, where smooth metal plates 
are at z = 0 and H. The periodic boundary condition is 
imposed along the x and y axes. We apply electric field 
under the fixed-potential condition, where the electric 
potential $(r) (defined away from the charge positions 
r ^ Vj) satisfies the metallic boundary condition, 

$ = 0 (z = 0), $ = -A$ = -E^H (z = H), (2) 

where A4> is the applied potential difference and E^, = 
A^/H the applied electric field. 

To realize the condition (2), we introduce image 
charges. For each charge Qj at rj = {xj,yj,Zj) in 
the cell, we consider images with the same charge qj 
at (xj,yj,Zj — 2Hmz) (rriz = ±1,---) and those with 
the opposite charge —Qj at (xj^yj, —Zj — 2Hmz) = 
0, ±1, • • •) outside the cell. Then, $(r) is expressed as 




m 
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\r - rj 


h\ 


■h\ 


Ea,^f 


( 3 ) 


where r ^ rj, fj = {xj,yj,—Zj), and the sum is over 
h = (Lmx, Lmy, 2 Hmz) with mx,my, and rriz being in¬ 
tegers. The sum over rrix and rriy arises from the lateral 
periodicity. Due to the sum over ruz, the first term in 
Eq.(3) vanishes at z = 0 and H, leading to 1^0) ■ The 
electrostatic energy Um is written as the surcPSlM]^ 





\rij +h\ 


<liQj 


-E^Mz, (4) 


where rij = - Vj, and = Vi - rj. In J^'ij^ 

exclude the self term with j = i ioi h = 0. The Mz is 
the z component of the total polarization Af, where 


M = (5) 

k i 


We can apply the 3D Ewald method due to the summa¬ 
tion over h fully accounting for the image charge^^^^^^. 

If the charge positions are shifted infinitesimally by 
dvi and Ea by dEa, is changed by 


dUn, = - ^ QiEi ■ dvi - MzdEz,, ( 6 ) 

i 


where Ei is the microscopic electric field acting on charge 
i = {k,a). Note that Ei consists of the contributions 
from the other molecules (I yf fc) for rigid water models. 
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The image charges are introduced as a mathematical 
convenience. The real charges are those in the cell and 
the surface charges on the metal walls. We write the 
surface charge densities as 170 ( 2 ;, y) at z = 0 and anix, y) 
at z = H. In our model, there are no charges in the 
thin regions 0 < z <1 A and 0 < H — z <l k due to 
the repulsive force from the walls, so we can set z = 0 
and H in Eq.(3) to find Eq.(2) (see a comment below 
Eq.(12)). The surface charge densities are then equal 
to Ez(x, y,0)/47r at z = 0 and —i?^(x, y, iJ)/47r at z = 
H, where Ez = —d^/dz is the electric field along the 
z axis. Their lateral averages a\ = f dxdy a\{x,y)/L^ 
(A = 0, H) are related to Mz 

do =-dH = E^/i-K + Mz/V. (7) 

The potential from the surface charges is divided into 
two parts as —Airdoz + (jJsi'f’), where 4 >s arises from the 
deviations Acto = ctq — (Tq and Aan = <Jh — dn- The 
microscopic field Ei in Eq.(6) is then divided as 

Ei = Ef + ATTdoBz + Es(vi), ( 8 ) 

where is the unit vector along the z axis. The Ef 
arises from the charges in the cell as 

^^ = EE' qjg{r^-rj+Lm^_), (9) 

3 

where mx = {mx,rny,0) and g{r) = r~^r. The ampli¬ 
tude of Ef is typically of order 15 V/nm and is very large 
for water (see Sec.III). The third term Es{r) = 
decays to zero away from the wallJ^, but it grows as 
the charges approach the walls (the image effect). If we 
neglect Eg retaining AirdoSz in Eg.(8), our method coin¬ 
cides with that of Yeh and Berkowitz^. 

The total potential U consists of three parts as 

U = Um + ^ y^»Lj(|2’fcO - 2-fol) 

^ k=jie 

T 'y ( [^w(-^fco) T ^ko')\i (i^i^) 

k 

where mlj is the pair potential and Uw is the wall poten¬ 
tial for the oxygen atoms expressed as 

UL 3 {r) = Ae[{a/rf'^ - {a/rf], (11) 

■Uw(z) = C' 9 (cr/z)® - Csia/zf. (12) 

We set e = 93.2fcB, a = 3.1589A, Cg = 27re/45, and 
C 3 = 15(79/2. The density n in liquid water is close to 
a~^ and cr characterizes the molecular length. Also due 
to the repulsive part of Uw, the molecules are depleted 
next to the walls, which leads to Eq.(6). 


B. Linear response and surface charge fluctuations 

The equilibrium distribution of our system is given by 
Z~^ exp(—"H/fcBY), where the Hamiltonian H = K. + U 


probability distribution 



FIG. 1: Gaussian distributions of (a = 

X, y, z) for A4> = 0 in (a) and 1.9 V in (b). Variance of SMz 
depends on A4> and is smaller than those of M^. and My. 


consists of the kinetic energy /C and the potential U in 
Eq.(lO) and Z is the partition function. The free energy 
is defined by F{T,Eg) = —kBTlnZ. From Eq.(4) H is 
expressed as H = Ho — E^^Mz, where Ho is the Hamil¬ 
tonian without applied field. Using this form, we may 
develop the linear response theory for small Ea. In our 
theory, the symbol (• • •) denotes the average over this 
distribution. However, in our simulation results, (• • •) 
denotes the time average over 6 ns. 

As in spin systems, the average polarization (Mz) and 
the polarization variance can be expressed a^^ 

= {i6Mzr)=kBT^^, (13) 


where SMz = Mz — {Mz) and the derivatives are per¬ 
formed at fixed T. For any physical quantity A (inde¬ 
pendent of Ea), the thermal average (A) is a function of 
T and Eg, and its derivative with respect to Eg, reads 


d 


{A) 


^(A<5M,). 


(14) 


As Ea —t 0, (A) is expanded with respect to Eg,, leading 
to the linear response expression. 


(A) = (A)o + {AMz)oEg,/kBT + ■■■, (15) 

where (• • •)o is the equilibrium average with Eg, = 0, so we 
have {Mz)o = 0. Using Eq.(7) we introduce the effective 
dielectric constant eeff of a film by 

Eeff = 47r(do)/Ua = 1 + 4^(M,)/UUa, (16) 

where V = L^H is the film volume. As Eg, — )■ 0, we find 

lim £eff = 1 + ^^{MDojVkBT. (17) 

In our previous papei^^, SeS from Eq.(16) in the linear 
regime was 21.4, which was close to that from the right 
hand side of Eq.(17). 

In Fig.l, we show that Mz, M^, and My obey Gaussian 
distributions. Here, {{5Mz)‘^)o/VkBT is 1.6 and 1.2 for 
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A$ = 0 and 1.9 V, respectively, while {M^)o/VkBT = 
{My)o/VkBT = 4.4 both for A$ = 0 and 1.9 V. Here, if 
A$ = 1.9 V, the applied field is Eg, = 0.42 V/nm. Thus, 
the fluctuations of are considerably suppressed than 
those of Mx and My. 

The surface charge densities consist of the deviation 
and the mean value. That is, we write erg = dg + 
A(Tg(a;, y) at z = 0. The mean dg gives rise to a ho¬ 
mogeneous electric field dTrdge^ in the cell as in Eg.(8). 
From Eq.(7), tfg is a fluctuating variable with variance^ 






kBT d{M,) 

E2 dEg, ’ 


(18) 


where Jdg = dg — (dg) is the deviation from the thermal 
average (dg). Thus, {{Sa'o)'^)L'^ kBT{eeff — l)/47riJ as 
Eg, —>• 0, which decreases as H~^ for large H. See Ap¬ 
pendix B for more discussions. In contrast, the deviation 
A(Tg(a:, y) i s mostly determined by the molecules close to 
the walP3 and the electric field Es{r) in Eq.(8) is no¬ 
ticeable only close to the wall. Recentl y, L immer et al. 
discussed the surface charge fluctuations^Sl. 


In our case, Stern layers are formed even for pure water 
with thickness about 5 A, where the molecules are under 
strong influence of the walls. In experiments, the Stern 
layers have been obser ved on ionizable dielectric walls in 
aqueous solution^^^^^, where the layers can be influenced 
by ion adsorption and/or ionization on the surface. 

Outs ide th e layers, a homogeneously polarized state is 
realizecPSEsHU] thickness H — 2d, where P{z) = Pg 
and —d'^{z)/dz = Pg with 

47r(CT-g) = Pg -I- 47rPg = ePg, (21) 

where e is the bulk dielectric constant. In Fig.2(a), we 
show Pb/u^o vs. A$. Here, n is the average density in 
the bulk (n = 1.071 /(t^ for A$ = 0 and n = 1.058/cr^ 
for A$ = 1.9 V). We can see that Pg increases linearly 
for A$ < 4 V but it tends to saturate for larger A4>. 
In our previous papei!^, we determined e directly from 
e = 1 -I- 47rPg/Pg to obtain e ^ 60 for A$ < 4 V, but 
this method is not accurate for very small Pg and Pg. 


2. Surface electric length 


C. Equilibrium states under applied field 

In equilibrium under applied field. Stern layers with a 
microscopic thickness d appear near the walls, where a 
potential drop occurs. However, dielectric response de¬ 
pends on another longer length in highly polar fluids. 
In this paper, the cell width Ed is much longer than d but 
can be shorter than 

1. Stern layers and bulk polarization 

We introduce a microscopic polarization density p{r) 
in terms of the charge positions Vj for polar molecules. 
It is related to the microscopic charge density p{r) bjEZl 


In terms of a surface electric length the sum of 

the potential drops in the top and bottom Stern layers 
is given by A$.^w/(^w + P), so the potential decreases 
by A$P/(£w -h P) in the bulk. As a result, the bulk 
quantities e and Pg are related to the film quantities £eff 
and Eg, by 

e/eeff = Pa/Pg = 1 + ^w/P- (22) 

In Fig.2(b), we plot £w vs. A4>, which is calculated from 
Eq.(26) below (so its accuracy is not good for very small 
A<I>). Her e, £w 9 nm for A<I> < 4 V. In the previous 
simulation^SEMI]^ ^ 2P. In our case, as Eg, —)• 0, we 
have iw/H = 1.86, e/ses = Pa/Pb — 2.86, and e ~ 60. 

From Eqs.(17) and (22), the variance of at Pa = 0 
is written as 


.p = p = Y^qj5{r-rj). (19) 

3 


{M!)o = 


VkBT f \ 

1 -h C/P W 47rP ) 


VkBTx 

1 + C/P’ 


(23) 


See Eq.(A4) in Appendix A for the expression of p. 
The total polarization (5) is given by the integral dSd = 
f drpir) in the cell. In our geometry, the average polar¬ 
ization P{z) = {pz{r)) along the z axis depends only on 
z. Using this P(z) and the average surface charge density 
(dg) we define the average Poisson potential 'k(z) by 

4'(z) = —47r((tg)z-I-47r f dz'P{z'), (20) 

Jo 

Here, 'I'(O) = 0, while p (H) = —HEg, follows from 
Eq.(7). The other group^^sElIII] determined 'k(z) from 
d^dt/dz^ = —47r(p)(z) without introducing p. For gen¬ 
eral geometries without ions, the average potential 4'(r) 
is obtained from = 47rV • (p(r)) with boundary con¬ 
ditions. 


where we can omit £w/47rP(^ x) for water. Hereafter, 

X = (e - l)/47r (24) 

is the bulk susceptibility. On the other hand, and My 
are decoupled from dg (see Eqs.(33) and (37)). Thus, 

(M2)o = (M2)o = VkBTx- (25) 

We express C in terms of P{z). From Eqs.(20) and 
(22) we obtain {Eg, — E\J)H = Att dz[Pg — P{z)], so 

C = (e - 1) f dz[l - P(z)/Pb], (26) 

Jo 

where the integrand 1 — P{z)/P\y is nonvanishing only 
around the Stern layers. Hence, is independent of P 
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FIG. 2: (a) Normalized bulk polarization density Ph/nno 

vs. AtO, where n/io is the maximum polarization, (b) Surface 
electric length £„ vs. Aik. Displayed also are profiles of 1 — 
P{z)/Ph at A<1? = 0.19 and 1.9 V near the bottom in (c) and 
the top in (d). Its integral is ^w/(e — 1) from Eq.(26). Due 
to depletion, Pb — 0 for 2 < 1.5 A and H — z < 1.5 A. For 
larger separation, variations for A"!* = 0.19 V are amplihed 
than those for Aik = 1.9 V because of smaller Pb. 

for ^ d. It is worth noting that the integration of 
1 — P{z)/Ph is analogous to the calculation of a Gibbs 
dividing surface between two coexisting phaseJ^. 

In Fig. 2 (c) and (d), we display 1 — P(z)/Pb for A$ = 
0.19 and 1.9 V near the bottom and top. Its integral is of 
order 1 A mainly due to the depletion at the walls, but 
£w is about 9 nm due to large e — 1. For A<i> = 0.19 V, 
we can see a maximum at 2 ^ 2 A in (c) and a minimum 
at — 2 ~ 2 A in (d), which originate from preference 
of the H-down configurations due to the image effeclP^. 
However, these extrema disappear for A$ = 1.9 V due to 
the influence of the increased surface charge density given 
by (d’o)/e = 0.44/nm^. Oscillatory behavior follows for 
larger separations. The curves ofA$ = 0.19V and 1.9 V 
are very different, but their integrals are nearly the same 
as in Fig.2(b). 

3. Surface capacitance 

The potential change in the Stern layers can be calcu¬ 
lated from the Poisson potential 'I'(z) in Eq.(20) as 

(A$)s = (do)/C+ + $- (2-O), 

= -(dff)/C_ - chS'o (27) 

where {o' h) = ~(dn) and C+ and C_ are the surface 
capacitance^^SmEMl] Due to the molecular anisotropy 


of water, there appears an intrinsic potential change at 
zero applied field given by <i>oQ = —47r Af dzP(z ), which 
sensitively depends on the wall potentiaP^IIIIlI]. In fact, 
was 0.2 V in Willard et al.’s simulatiorl^, but it was 
—0.09 V for weaker adsorption in our simulatiorP^. 

The sum of these potential drops gives the total bound¬ 
ary drop in the form ag/C with 

C = + CZ^)-\ (28) 

where does not appear. In terms of this C, the length 
is expressed as 

= e/AnC. (29) 

It is worth noting that the surface capacitance of a di¬ 
electric film on a metal wall is given by Cd = EdlPnid^ 
where Ed is its dielectric constant and Id is its thickness. 


D. Polarization correlations 

1. Dipolar and homogeneous correlations 

We consider the polarization correlation function, 

Ga/3(ri,r2) = {5pa{ri)Spp{r2)), (30) 

where a,/3 = x,y,z and Spa = Pa — {Pa)- In this sub¬ 
section, En is small in the linear response regime, so 
(•••) = (•• -A in Eq.(30) and ((<5M„)2) ^ 

FelderhoP presented a continuum theory of the polar¬ 
ization fluctuations for infinite, uniform systems without 
electrodes. Eor isotropic polar systems, he found short- 
ranged and dipolar correlations as 

Gap{ri, r2)/kBT = xSa/36{r) + /e)V i3r~^, (31) 

where r = r 2 — Vi and Vq, = djdxa- Here, = 

2 iXaXj 3 jr^ — 5ap jr^^ which diverges as r —>■ 0 and is mean¬ 
ingful only for r > a. The delta function in Eq.(31) 
should also be treated as a normalized shape function 
with a width longer than a. We may conveniently remove 
this divergence of the dipolar term by replacing r~^ by its 
long-range part, written as ifi jr). I n the Ewald method, 
use has been made of the fornP^Hil^ 

ipe(r) = erf ('yr)/r, (32) 

where 7 ^ and erf(M) = {2/^/tt) due~'^^ is the 

error function. Here, ipi = r~^ for r ^ 7 “^ and 
■0^ = ( 27 /y^)(l — ■ •) for r <C 7 “^. Notice that 

the inner product Gq,q, = {Sp{ri) ■ 6p{r2)) has no 
long-range dipolar correlation. The dipolar correlation 
itself follows in the continuum theory (from Eq.(Bl) in 
Appendix B), as in the case of dipolar ferromagnetP^. 

In our theory, Eq.(7) indicates that the fluctuations of 
the mean surface charge density do produce polarization 
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fluctuations homogeneous in the bulk along the z axis 
(see Appendix B). We propose the following form, 

Gaii{ri,r 2 )/kYiT = x5afi5{r) + {x^/e)V p'ipi{r) 

-\-5az5pzho/V, (33) 

where ri and r 2 = ri +r are in the bulk and ho is a con¬ 
stant to be determined in Eq.(39). We have replaced r~^ 
by tpi and neglected deformation of the dipolar correla¬ 
tion due to the metal walls. If we pick up two molecules 
k and i in the bulk, their correlation at large separation 
(t) tends to a constant as 

(/Xfc • H() = ■ ne) = k-QTho/v?V, (34) 

where the dipolar correlation vanishes. 

Setting n = ro -I- r' and r 2 = ro -I- r in Eq.(33), we 
integrate Gaa = Gaai'i'o + r', Tq -|- r) {a = x, y, z) with 
respect to r = (x, y, z) within a spheroid expressed by 

ix^ + y^)/Rl+zyRl<l, (35) 

where the radii R± and i?|| much exceed 7 “^ and r' is 
also within the spheroid (away from its surface). The two 
points Tq + r' and + r are in the bulk region. These 
integrations can safely be performed because the dipolar 
terms — r'|) are finite at \r — r’\ = 0. The 

resultant integrals are independent of r' aJ^ 

/ drG,,lkBT = x - X^N, + ho^, (36) 

J spheroid ^ * 

r Ott 

/ drG^^lk^T = x -x'(l - N,), (37) 

d spheroid ^ 

where v = 47 ri?||i?i ^/3 is the spheroid volume. The in¬ 
tegral of Gyy is equal to that of Gxx- The N^, is the 
depolarization factor along the z axis determined bj 0 ^ 

77 = i?ll/i?_L. (38) 

The Nz tends to 1 for 77 ^ 1, 1/3 for 77 = 1, and 0 
for 77 ^ 1. The last term in Eq.(36) arises from the 
homogeneous term in Eq.(33), so it is proportional to v. 

We also integrate G^^(ri,r 2 ) in Eq.(33) with respect 
to ri and r 2 in the whole bulk region, where we can use 
Eq.(36) in the pancake limit = 1 under the lateral 
periodic boundary conditiorP^. In Appendix B, we shall 
also see that the polarization sum in the Stern layers is 
much smaller than that in the bulk region. Thus, we 
have {{6Mz)‘^)/k^T = I4(x/e + ho)- From Eq.(23) we 
now determine ho as 

ho = x/(l + £w/H) - x/e = (Eeff - l)x/e, (39) 

which considerably depends on H for H < but tends 
to 47rx^/e for H Another derivation of Eq.(39) 

will be presented in Appendix B. For A/, = 1 in Eq.(36), 
the first short-range term and the second dipolar term 
almost cancel for £ 1 , but the third homogeneous term 

increases with increasing v leading to Eq.(23) for v = V 
(see Fig. 14). We also obtain Eq.(25) from Eq.(37) with 
Nz = 1 , where the dipolar part does not contribute. 


A(I)=0 



FIG. 3: Gaussian distributions of M%/{vk-BT)^G (q, = 

X, y, z) in a spherical cavity, whose radins R is (a) 0.95 nm and 
(b) 1.58 nm at A4> = 0. In (a), there is almost no difference 
in the three directions with variance ((M® )^) = 2.9vkBT. In 
(b), it is decreased to 2.3vkBT for a = z and is increased to 
and 3.2vkBT for a = x,y due to the boundary effect. 

8. Kirkwood-Frohlich formula 

To reproduce Kirkwood-Frohlich formula, we define 
the polarization integrals SM^ = J^^j^drSpa{ro + r) 
within a sphere with radius R. Then, we find 

{{6Ml)^)/vk^T = x( 2£ + l)/3e + dozhovIV. (40) 

Without the second term, the above relation is well- 
known in the literatur^i^. For a sphere in infinite sys¬ 
tems, the factor (2£-|-1)/3£ in the first term was shown to 
arise from the reaction field produced by the molecules 
in the exterioi!^. In Appendix C, we will present an¬ 
other simple derivation of Eq.(40) in the limit v/V —)■ 0. 
In Fig.3, obey Gaussian distributions. For R = 

0.95 nm, ((M®)^) = 2.9vkBT for the three directions 
consistently with Eq.(40). However, the boundary effect 
appears for R = 1.58 nm, where ((M|)^) = 2.37 ;/cbT and 
((M|)2) = ((M^)2) 7^ 3.2vkBT. 

The left hand side of Eq.(40) has been written as 
ny^gY^/kBT for infinite systems (v/V —>■ 0) in terms of 
the correlation factor ^k, which leads to the Kirkwood- 
Frohlich formula for nonpolarizable molecule^^^, 

X = {nylgK/kBT)e/{2e + 1). (41) 

The Onsager formulcP follows for (/k = 1. The factor 
(7 k is defined in terms of the dipole correlation between 
molecules k and i in the bulk apl 

9K = ^ + '^{nk-ne), (42) 

i 

where the first term 1 arises from the self term (i = k) 
and the sum in the second term is over ^(^ k) at fixed 
k with molecule £ being in a sphere with volume v <^V. 
Here, the dipolar correlation vanishes and the nearest 
nei ghbo r molecules give a dominant contribution in the 
sunJlHl. This assures that the formula (41) holds in the 
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thermodynamic limit. This has been confirmed by many 
authorpmHU. 

In our case, we obtain = 2.14 from Eq.(42) at 
A$ = 0 (see Fig.15(b) in Appendix D), which yields 
X = 4.68 and e = 59.8 from Eq.(41). Now, in addition to 
this derivation, we have determined e from the dielectric 
relation (21)23 and from the data of (M|)o and (see 
the sentence below Eq.(22)), all yielding e = 60 in the 
linear regime. A close value of e was obtained by Abas- 
cal and VegaP^for the TIP4P/2005 model under the 3D 
periodic boundary condition. 

It is worth noting that Harris and Aldei® added the 
term najn{e + 2)/3 in the right hand side of Eq.(41) for 
polarizable fluids, where am is the molecular polarizabil¬ 
ity. Their expression reduces to the Clausius-Mossotti 
formula 3x/(e + 2) = na^ as fiQ —>■ 0. It was used to 
analyze a wide range of data on e for wateJ^. 


III. NUMERICAL RESULTS 

In this section, we further present results of our molec¬ 
ular dynamics simulation on the fluctuations of the local 
fields and the dipoles. 


A. Fluctuations of local field 

Electric field fluctuations have been nume rically exam¬ 
ined in water and electrolytes (for A4> = o) 27 | 33| t35]^ 
cause they play crucial roles in dissociation reactions and 
vibrational spectroscopic response. Here, we examine the 
thermal fluctuations of the local electric field which 
is defined by Eq.(A2) in Appendix A. 


1. Joint distribution 

For each molecule fc, the direction of is written as 
efe = \E'n-^E'r. (43) 



FIG. 4: (a) Joint distribution P]oc{E,cos8) of the local 

electric field for E = and cos6» = Eiiy\E^^^\ in 

Eq.(44), where A<F = 1.9 V. (b) Conditional probability 
distribution Pioc{E, cos 9)/P{^{E) and (c) its logarithm for 
E — 13.0,16.8, and 20.6 V/nm, where P{!^{E) is the distribu¬ 
tion for E — in Eq.(46). This conditional distribution 

is of the exponential form oc exp(Rcos0). (d) B depends on 
E linearly. 


In terms of its z component Czk, the angle of the local 
field with respect to the z axis is given by cos“^(ezfe). 
We consider the joint distribution for and Czk- 

Pioc(E, cos e) = ((5(|£;L°"| - E)Siezk - cos 0))b, (44) 

which is related to the 3D local-field distribution as 

P®(£;) = {SiE^°^-E))^ 

= {27rE^)-^Pioc{E,cos0). (45) 

Hereafter, (• • •)b = (X^feebuik-^b"^ ’'') denotes the aver¬ 
age over the molecules in the region 0.3H < ZkG < 0-7H 
and over a time interval of 6 ns, where zj-g is the z com¬ 
ponent of the center of mass of molecule k and Ab is the 


number of the molecules in this region. The distribution 
of and that of Czk are separately written as 

PitiE) = dsPUE, s) = (5(|£;'r I - E))^, (46) 

POO 

Pi” (s) = / dPPioc(P, s) = {S(ezk - s))b, (47) 

Jo 

where we write s = cos0. 

In Fig.4, we display (a) Pioc{E,cos9) and (b) the con¬ 
ditional probability distribution Pioc{E, cos 9) / P^{E) in 
the E-cos 9 plane at A4) = 1.9 V. The former is peaked at 
an intermediate Eq independent of cos 9, while the latter 
depends on E very weakly in the displayed range of E. 

















AO=0 



FIG. 5: (a) Distribution Pioc(i^) for E = in Eq.(46) at 
Ail? = 0 (red) and Gaussian distribution P^{E) in Eq.(50) 
(blue). The two curves nearly coincide, (b) Ratio of this dis¬ 
tribution to the Gaussian one (= P^{E)/(E)) at A<E> = 0. 
Shown also is the ratio of the distribution at A<1? = 1.9 V 
to that at A<1? = 0 (= Pioc(P)a«=i. 9 v/Pioc(P)a«=o) (in¬ 
set). They are very close to 1 in the displayed range 9 V/nm 
< E < 25 V/nm. 


In (c), the latter is nearly of the exponential form, so 
P (F rn-P) ~ pab/ p^ exp[P(P)cOS0] 

Pioc(P) COS t'j — Piq^[E) Zi^BiE')') ’ 


where Z{x) = 2 sinhz/x is the normalization factor. The 
coefficient B = B{E) weakly depends on E as 


B[E) ^ B{Eq) + 0.46(P/Po - 1), (49) 

where Eq gives a maximum of Pi^{E) (see Eq.(50)). 


(b) snapshot including 

(a) H^O and 0->H bonding second nearest molecules 



E [V/nm] E [V/nm] 


FIG. 6: (a) Typical configuration of 4 hydrogen-bonded 

molecules (first-nearest ones) around a reference molecule at 
the center. Hydrogen bonds (black bars) are from oxygen (hy¬ 
drogen) of the reference one to hydrogen (oxygen) of the other 
ones, written as O—?H (H—?0). (b) Snapshot of hydrogen- 
bonded molecules composed of 4 first-nearest ones and 11 
second-nearest ones, (c) Distribution Pit{E,K) in Eq.(53) 
for local-field contributions from group K: (i) O—?H, (ii) 
H—?0, (iii) first(0—?H+H—?0), and (iv) first-l-second ones. 
They are compared with P{^{E). (d) Distribution P{"^^{E, K) 
in Eq.(54) for local-field contributions from group P', where 
the contributions from group K are excluded. 


2. Large local-field amplitude due to hydrogen bonding 


In Fig.5, P^^{E) is almost independent of A$ in the 
linear response regime and is very close to a shifted Gaus¬ 
sian distribution expressed as 


P^^{E) = 




(E-Eq)^ - 

2so 


(50) 


whose mean and standard deviation are given by 

Eo = 16.8V/nm = 1.2e/cr‘^ = Sl/csT/zio, (51) 
= 3.8V/nm = 0.23Po- (52) 

Sellner et calculated P{^{E) at O and H sites for 
A$ = 0 to find similar behavior. 

In water, it is natural that a large fraction of the lo¬ 
cal field is produced by the hydrogen-bonded molecules 
around each molecule. Reischl et al^ analyzed this as¬ 
pect at the centers of the molecules’ OH bonds. Note 
that there are a number of definitions of the hydrogen 

bondPHm 

As in our previous papei!^, we treat two 
molecules to be hydrogen-bonded if one of the intermolec- 
ular OH distances is shorter than 2.4 A and the angle be¬ 
tween the 00 vector and one of their intramolecular OH 


bonds is smaller than 7r/6. A similar definition was used 
by ZielkiewicJ^. The average number of the intermolec- 
ular hydrogen bonds is then 3.6 per molecule in the bulk. 
For this HO-distance definition, the hydrogen bonds for 
each molecule k extend either from its protons (H—>-0) 
or from its oxygen atom (O—?>H). This classification is 
convenient for the present problem. 

In Fig.6, we display a typical configuration of 4 
hydrogen-bonded molecules (first nearest ones) in (a) and 
a snapshot of 15 ones (first and second nearest ones) in 
(b) around a reference molecule at the center. The second 
nearest ones are those hydrogen-bonded to the first near¬ 
est one^^. These surrounding ones give rise which 

is nearly parallel to (see Fig.9 below). In Fig.6(c), we 
write the distributions, 

- A))b. (53) 

For each reference k, are the local-field contri¬ 

butions from the other molecules in group K, where K 
represents (i) O—/-H, (ii) H—>-0, (iii) first(0—/-H-l-H—>-0), 
and (iv) first-l-second ones. To P^{E), the distribu¬ 
tion from the first nearest ones is considerably close and 
that from the first-fsecond ones is very close. Also the 
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local-field contributions from the nearest H—>0 bonded 
molecules are most important. This is because is 
rather close to +.£'fcH 2 ) (see Eqs.(A2) and (A3) 

in Appendix A). Furthermore, in Fig.6(d), we display the 
local-field contributions from group K which consists of 
the molecules not belonging to group K. Their distribu¬ 
tions are written as 

Pif,(F, k) = {5{\kr - E^m - E)k. (54) 

Naturally, the distributions excluding the first nearest 
neighbors and the first-l-second ones have smaller mean 
values and variances. 

3. Orientation distribution of local field 

In Fig.7(a), the angle distribution P[“(cos0) is excel¬ 
lently fitted to the exponential form. 



FIG. 7: Distribution P°J^{cos6) for the local field orienta¬ 
tion cos 9 = Czk (left) and that Pd (cos 0) for the dipole ori¬ 
entation cos9 = Uzk (right) for A>I> = 1.9 and 3.8 V. They 
nearly coincide and their logarithms are linear in cos 6 (inset), 
leading to the exponential forms in Eqs.(55) and (59) with 
Pioc = 0.85(1.71) and Bd = 0.84(1.68) for A4> = 1.9 (3.8) V. 


Pi”(cos6») ^ exp(PiocCos6')/Z(Pioc). (55) 

The coefficient B\ac is related to B{E) in Fq.(48) by 

Bloc = B{Eo). (56) 

The above form also follows from integration of Fq.(48) 
with respect to E with the aid of Fqs.(49) and (50), 
where we use the sharpness of Pif^{E) or the inequal¬ 
ity so/Eq ^ 1. See Fig.8(a) for Pioc vs. A$, where 
the linear growth (oc Ad)) occurs for Ad* < 4 V as in 
Fig.2(a). From Fqs.(48), (50), and (55) the bulk average 
of the z component of the local field is written as 

Floe - (Fife")b = Fo(e,fc)b = Fo/:(Fioc), (57) 

where C{x) = coth(a;) — 1/a; is the Langevin function. In 
the linear regime Pioc 1, we have Fioc = FqFioc/S. 

We may also consider the distribution of the l ocal fi eld 
in one direction, which was calculated previouslji^^HHUl 
It is defined by PiJ^{Ez) = {^{E^zk ~ Ez))h along the z 
axis. From Fqs.(45), (48), (50), and (55), we obtain 

p.D,p,..expW„.E./Bo) r^pP£,(E') 



FIG. 8: (a) Coefficients Pioc and Pd in the exponential 

orientation distributions in Eqs.(55) and (59) vs. Ad?, (b) 
Pp/Pb vs. Ad>, where Pp is the effective field proportional to 
Pd in Eq.(61) and Pb is the bulk field in Eq.(21). 

This emonential form was assumed in the classical 
theories!)^. Recently, it was numerically obtained for 
water by He et alW^. Furthermore, in Fig.7, Pd(cos0) is 
very close to P[°/.(cos0) for the local field orientation in 
Fq.(55). In Fig.8(a), the difference between Pioc and Pd 
is indeed very small. In terms of Pd, the average bulk 
polarization Pb in Fq.(21) is expressed as 


For A$ = 0, this distributio n exh ibits a plateau in the 
range \Ez\ < Eq from Fq.(50)P3^. For Ad? = 1.9 V, its 
profile was calculated in our previous papei!^. 

B. Dipole fluctuations 

1. Classical exponential form 

For the polarization direction rik = we consider 

the distribution of its 0 component Uzk = cos 9k- In 
Fig.7(b), it is again nearly of the exponential form, 

Pd(cos6») = {6{nzk - cos9))b 

= exp(Pd cos0)/F(Pd). 


Pb = nfj.o{nzk) = nfj.oC{Bd). (60) 

Here, we introduce the effective electric field Pp 

Pp = kBTBd/fio = (fcBT//ro)F“^(Pb/n/ro), (61) 

where C~^{x) is the inverse function of L{x). In the 
linear response regime, we obtairf^ 

Pp = SkBTPh/nfil = [ieg^/ (2e -b l)]Pb- (62) 

In Fig.8(b), we plot Pp/Pb using the data in Fig.2(a), 
which is about 3 in the linear regime but decreases con¬ 
siderably in the nonlinear regime. 

Since Pioc = Pd, Fqs.(57) and (60) give 

Pioc/Pb = Eo/nfio = 7.0. 


(59) 


(63) 
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A<D=0 



FIG. 9: (a) Distribution Pm (cos (/j) for cos(fi = rik ■ e.k in 

Eq.(66) and P™(cos(p) for cosi/p = fik ■ e.^ in Eq.(68) at 
Ail? = 0, where fik is the time-averaged polarization direction 
in Eq.(67). Their logarithms are also shown (inset). Then, 
Pm(cos(/?) oc exp(19.6 cos 99), but P™(cos9p) does not exhibit 
the exponential form for small ip. (b) Pid(cos ip) sin and 
Pi){"(cos (p) sin vs. p> at A<1? = 0. 


Here, we introduce the Lorentz factor 710 c by 

E\oc = Ex, + dTryiocPb- (64) 

From Eq.(63), 710 c can be expressed in terms of Eq as 

7ioc = {Eo/fion - l/x)/ 47 r, (65) 

which gives 710 c — 0.56 even in the nonlinear regime. 
In our previous papeil^, we directly calculating E\oc = 
{E^°k)b to obtain 710 c = 0.58 for any A$. Here, if we set 
Pb = n{am+fJ-l/SkBT)Ep wit hPp = Eioc and 710 c = 1/3, 
we obtain the Debye formula!^ 3 %/{e + 2) = n{am + 
/ig/3A:BP)) where am is the molecular polarizability. For 
liquid water, however, we find Pioc/^'p — k-oEo/^kBT = 
10 from Eqs.(62) and (63) in the linear regime. 


2. Nearly parallel nt and et 


AO=0 



FIG. 10: Distributions for vector between the two protons 
in Eq.(69) and its conjugate force P), in Eq.(70). (a) 
in Eq.(71) for the amplitude at A$ = 0, which is 

Gaussian form with mean 10.8 V/nm and standard deviation 
3.1 V/nm. (b) P/’''(cosx) in Eq.(72) for the angle between 
^ and Fk at A<1? = 0, where these two directions are nearly 
parallel. Shown also is its logarithm (inset), indicating the 
form oc exp{B^ cos 9) with = 16.6. 

changes on timescales of order 5 ps. Thus, we define the 
temporally smoothed polarization direction, 

/ At /2 

dTnk{t + T), (67) 

-At/2 

where we set At = 0.03 ps and Afk is the normalization 
factor ensuring \fik\ = 1 with |(A/’|)b|^^^ — O.OAt. The 
angle (pk = cos~^{fik ■ e^) for fik should be smaller than 
ipk- In Fig.9(a), we plot its angle distribution, 

PiX(cos(/?) = ((5(cos(/? - fik ■ efc))b, (68) 

which is surely sharper than Pid(cos:/?) in Eq.( 66 ) with 
(cos(^fc)b = 0.97 and ((cos(^fc)^)b = 0.94. In Fig.9(b), 
the curve of Pm™(cos ip) sin ip vs. (p is more sharply peaked 
for small tp than Pm (cos (p) sin (p. 


The result Pd (cos 0) = Pi“(cos0) in applied field in 
Fig.7 indicates that the dipole direction and the local 
field direction should be nearly parallel both without 
and with applied field, despite large thermal fluctuations 
at P = 298 K. As discussed in Sec.I, this is because of 
large 30 in Eq.(51). We thus consider the 

angle (pk = cos“^(nfc • Sk) between the two directions. In 
Fig.9(a), its distribution at A$ = 0 behaves as 

Pid(cos(/?) = ((5(cos(/? - Uk ■ ek))b 

= exp(PidCOS(/))/Z(Pid), (66) 

with Pid = 19.6. We obtain {cos(pk)h = 0.95 and 
((cos(/)fc)^)b = 0.90, so y(Pfc)b ~ 0.3 ~ O.Itt. In 
Fig.9(b), we also plot Pm (cos:/?) sin </) vs. p, which is 
approximated as PM 7 ?exp(—P mP^/ 2) for small ip. This 
distribution is insensitive to A$ (not shown here). 

Moreover, as will be shown below, rik (t) contains a fast 
librational part on timescales of order 0.01 ps, while ek{t) 


C. Relative vector between two protons 

For each molecule k, we may also consider the relative 
vector between the two protons, 

Ifc = '^feHl —(69) 

whose length is fixed at = ohh = 1.514 A. From 
Eq.(Al) in Appendix A, its conjugate force is given by 

E\ = 2 ®h(£^/cHi — EkB2)- (70) 

As in the dipole case, we first consider the distribution 
for the field amplitude iP^I/gn = — £lfcH 2 |/ 2 : 

pf(p-^) = (5(p^-ini/9H))b. (71) 

In Fig.10(a), P^'°{E'‘) is of the Gaussian form in Eq.(50). 
Its mean value and standard deviation are 10.8 V/nm 





















11 


and 3.1 V/nm, respectively. This form indicates strong 
asymmetry between the forces acting on the two protons 
due to the hydrogen bonding. For a large-angle rotation 
of the energy needed is of order aHH(|.F’5.|)b = 29kBT. 
Thus, and F]. should be nearly parallel in equilibrium, 
as in the case of and In Fig.10(b), for the angle 
Xk = cos“^(^j;, • F^/ohhI.F’^I) between these vectors, its 
distribution is again of the exponential form, 

P”(cosx) = ((5(cosx - cosxfe))b 

^ ^ cos x)lZ{B^), (72) 

where B^ = 16.6, so (cosxfc)b = 0.93 and ((cosxfe)^)b = 
0.87. The surely tends to be parallel to F\. 

D. Orientation time-correlation functions 

We examine the single-molecule time-correlation func¬ 
tions averaged over the molecules in the bulk. First, we 
consider those for the local field direction efc(t) and the 
dipole direction rikit)-. 

C’ioc(t) = (efc(t)-efc(0))b, (73) 

Cd(t) = (nfc(<) • nfc(0))b. (74) 

In Fig.11, we plot C\oc{t) in (a) and Cd{t) in (b), which 
are nicely fitted to the exponential form oc exp(—t/rh) 
with Th = 5.0 ps for t > 1 ps. This indicates that nk{t) 
moves with (t ), where e^, (t) should be gov erned by the 
hydrogen bond dynamics. In the literatur cj^^ l ^^ l ^^ l ^ -^l^, 
Cd{t) has been found to decay exponentially for t > 1 ps 
at r ~ 300 K, while it is known to exhibit a str etched- 
exponential relaxation in supercooled state^^^^. 

Furthermore, we notice that C'd(t) has a small local 
minimum at t = 0.03 ps. To understand this short-time 
behavior, we consider the perpendicular component of 
nfc(t) with respect to efc(t) defined by 

nkj_{t) = nk{t) - [nfc(t) • efc(t)]efe(t), (75) 

which satisfies nk^_{t) ■ ek{t) = 0. In Fig.11(c), we plot 
its single-molecule time-correlation function, 

C±(t) = (nfc_L(t) • nfej_(0))b, (76) 

where C'_l(0) = (|nfc_Lp)b — 0.093. We recognize that 
nk±{t) undergoes librational motions and C±{t) decays 
on a timescale of 0.01 ps with a minimum at t = 0.023 
ps. In liquid water, similar short-time behaviors with lo¬ 
cal minimum and maximum have been observed in vari¬ 
ous time-dependent quantities in simulations and exper¬ 
iments as manifestation of librational motion^2lHlI]_ 

To examine the long-time and short-time behaviors of 
C'd(t) separately, we also define the parallel component 
’^fcll(^) = [^fc(^) ■ efc(t)]efc(t). Then, we have the decom¬ 
position Cd{t) = C'ii(t) -I- C±{t) + C'cr(t) with 

CyW = («/c||(i) •«fc||(0))b, (77) 

C'cr(t) = {nk\\{t) ■ nfc_L(0))b -f {nk±{t) ■ nfe||(0))b, (78) 


A(1)=0 



t [ps] t [ps] 


FIG. 11: Single-molecule orientation time-correlation func¬ 
tions. (a) Cioc(t) for the local field direction ek{t) and (b) 
C'd(t) for the dipole direction rikit) (red lines), where ex¬ 
ponential functions Texp(—t/rh) are written with common 
Th = 5.0 ps (dashed lines), (c) C±{t) for the perpendicu¬ 
lar component nk±it) in Eq.(75), which exhibits an oscil¬ 
latory decay on a timescale of 0.01 ps. (d) Decomposition 
C'd(t) = C\\{t) + C±{t) + C'cr(t) using Eqs.(56)-(58) with ex¬ 
ponential htting for C\\ (t) (dashed line). Short-time minimum 
of Cd{t) in (b) arises from that of C±{t) in (c). 

where C'||(0) = 0.907 and C'cr(O) = 0. In Fig.11(d), we 
plot these three time-correlation functions. The domi¬ 
nant term C\\{t) has no short-time minimum and decays 
exponentially for t > 1 ps. Thus, the short-time min¬ 
imum of Cd(t) arises from the rapid motions of nk±it). 
The cross contribution Ccr(t) is at most 0.08 decaying 
exponentially for t > 1 ps. 

The scenario in Fig.11 holds also for the relative vector 
between the two protons in Eq.(69) and its conju¬ 
gate force i^fc(t) in Eq.(70). In Eig.12(a), displayed are 

Cdt) = {Ut)-U0))/<H, (79) 

Cfdt) = {fkit)-fki0)), (80) 

where fdd = represents the force direc¬ 

tion. These functions decay exponentially as exp(—t/r^) 
with Tj = 5.8 ps, which is slightly longer than Th = 5.0 
ps in Eig.11. To detect librational motions, we present 
the time-correlation function, 

Q_L(i) = (4fe_L(i) • 4fe_L(0))/aHHj (81) 

for the vector =^kit)-[^kd)-fk{t)]fkd) perpen¬ 
dicular to /fc(t). In Eig.12(b), closely resembles 

C±{t) in Fig. 11(c), decaying on a timescale of 0.01 ps. 
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Acl>=0 



FIG. 12: Single-molecule orientation time-correlation func¬ 
tions related to relative vector between two protons, 

(a) Cf/^it) = {flit) ■ fl{0)) and ^(t) = ■ $fc(0))/a^H, 

where fl{t) = |F’5c(^)l These exhibit long-time ex¬ 

ponential decay oc exp(—t/r^) with common = 5.8 ps. (b) 
C^±(t) for vector ■ fl{t))fl{t)] perpendic¬ 

ular to flit), decaying rapidly as C±it) in Fig.11(c). 


A(1>=0 



FIG. 13: Detection of long-range dipolar correlations in 

{ro,ro+r)r^ for Atl? = 0, where the self part is excluded, 
(a) jk^T , where cos Q = zjr is fixed at ±1, ±l/\/3, and 

0. (b) Its value at r = 10 A, which is fitted to 0.4(3 cos^ 9 — 1) 
(broken line), (c) Combination X]q /3 XaXpr/k^T, whose 
limiting value is about 0.76. (d) Trace Glfffr^/k^T con¬ 
taining no long-range dipolar contribution. 


The short-time librational motions are appreciable for 
the protons but are small for the oxygen atoms, because 
of the large mass ratio. In fact, the mean square displace¬ 
ment of the protons (|rf;Hi(^ + At) — r/jHi(t)P)b exhibits 
a plateau-like behavior at 0.2 around At ^ 0.04 ps. 
This aspect will further be studied in future. 


AO=0 



FIG. 14: 7^2 and + lyy)/‘2 in Eq.(82), which are integrals 
of Gaa{ri, r 2 )/k^TL^b with respect to ri and r 2 in the plate 
region with thickness h (light blue region). Their behaviors 
for b ^ a are predicted in Eqs.(83) and (84). For b > a, 
Eq.(83) holds nicely, but Eq.(84) holds only approximately. 


E. Polarization correlations 

We calculate the polarization correlation function Gap 
in Eq.(30). There have been some attempt^^sHsU 
tect the long-range dipolar correlation in liquid water, 
but they have not been based on Felderhof’s expression 
(31). We aim to detect it unambiguously in agreement 
with his theory. Note that it vanishes in the radius- 
dependent correlation factor Gk{R) in Eq.(D6), which 
has been calculated by many authors. We also detect 
the homogeneous term in Eq.(33). 

To detect the dipolar correlation, we consider the 
pair part omitting the self part in GapirijVi + r) 
(see Appendix D). For simplicity, we use a point-dipole 
approximatiorP^, which is not accurate for small r but is 
allowable for large r. In Fig. 13, we calculate the prod- 
uct From Eqs.(31) or (33), it should behave as 

ix'^/e)i3xaXp/r'^ — 1) for r a with — 0.37 (from 
our data of ^k)- In (a), we plot /k^T for A$ = 0, 

where cos6> = zjr = ±l,±l/-\/3, and 0. In (b), the 
limiting value can indeed be fitted by 0.4(3 cos^ 6* — 1). 
In (c), the combination 'Yliap^^ap tends to 

a constant about 0.76, while its theoretical value is 
2x^/e = 0.74. In (d), the trace X)a is dis¬ 

played, which exhibits no long-range dipolar contribu¬ 
tion. If it is divided by jk^T and is integrated with 
respect to r in the range 0 < r < 7?, we obtain Gk(I?) — 1 
(see Eq.(D6)). On all the curves in Fig. 13, we can see 
peaks corresponding to nearest neighbor molecules. 

To detect the homogeneous term in G^z in Eq.(33), we 
integrate G^z. over a plate region parallel to the xy plane 
with thickness b, where the 2 coordinate is in the range 
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[{H — b)/2, [H + b)/2]. We integrate Gaa{'i'i,r 2 ) with 
respect to ri and r 2 in this region. The integrals are 
divided by k^^TL'^h to give 

Iaaib)= f dri [ dr 2 Gaa{ri,r 2 )lk^TL^b. (82) 

J plate J plate 

In the pancake limit in Eqs.(36) and (37) in the lateral 
periodic boundary condition {N^ = 1), we obtain 

hz = x/e + Kb/H, (83) 

IXX = lyy = X) (84) 

which are valid for 6 ^ cr. The first relation also follows 
from Eq.(BlO) in Appendix B. In Fig.14, the curve of Izz 
vs. b nicely satisfies Eq.(83) with ho = 1.59 for 5 A < 
6 < 40 A, but it exhibits a small peak at 6 = 41 A. The 
curve of {Ixx + Iyy)j‘2 tends to x at 6 = but exceeds 
X by 10% for 25A < b < 40A, where the latter behavior 
is consistent with that of and in Fig.3(b). Both 
Izz and {Ixx + Iyy)/‘^ decrease for 40A < b < H due to 
the depletion layers (see Sec.IIA). 

IV. SUMMARY AND REMARKS 

We have studied the thermal fluctuations of the lo¬ 
cal field and the dipole moments in liquid water in ap¬ 
plied electric field with analytic theory and molecular 
dynamics simulation. In the latter, we have used the 
TIP4P/2005 modeP^ and the 3D Ewald methocPlHIZl jn 
a system of 2400 molecules between metal walls at z = 0 
and H = 44.7 A. We have included the image charges 
and fixed the potential difference A$ = HEg,. In the 
following, we summarize our main results. 

(1) We have derived linear response expressions from 

Eq.(4) in Sec.IIB. In particular, the effective dielectric 
constant eeff of film systems is expressed in terms of 
the polarization variance {M‘^)q at zero applied field in 
Eq.(17), which is smaller than (Mj)o = {My)o by the 
factor {1 + Here, M is the total polarization. 

The is a surface electric length (~ 10 nmjp^. 

(2) In Sec.IID, we have presented a theory on the po¬ 
larization correlation function Gap{xi,r 2 ) in Eq.(30) in 
the presence of metal walls. Together with the dipolar 
terrrP, we have found a homogeneous correlation term 
stemming from the surface charge fluctuations. In Ap¬ 
pendix B, this global coupling has been examined for the 
quantities averaged in the xy plane. It is important for 
the fluctuations of the total polarization Mz- However, it 
is negligible in a small specimen region in the bul k, s o we 
have reproduced the Kirkwood-Frdhlich formulsPEl. In 
Appendix C, we have presented another derivation of this 
formula, where an increase of the electrostatic energy is 
calculated for the polarization fluctuations. In Sec.HIE, 
we have numerically derived the long-range dipolar term 
and the homogeneous term in Gap- 

(3) In Sec.HIA, we have examined the fluctuations of 
the local field in applied field. The distribution P\oc{^) 


for its amplitude has been found to be of a Gaussian 
form in Eq.(50) with a large mean Eq = 17 V/nm and a 
relatively small variance. This form is mainly due to the 
hydrogen-bonded molecules, as illustrated in Fig.6. The 
joint distribution P\oc{E, cos 9) for the amplitude and the 
direction of the local field is then equal to the product 
of T’ioc(A) and the exponential one oc exp[H(A) cos 0], 
where B{E) depends on E linearly as in Eq.(49). The 
orientation distribution Pj°(,(cosd) is of the exponential 
form oc exp(i?ioc cos 9) with B\oc = B{Eq). 

(4) In Sec.HIB, the distribution P^{cos9) for the po¬ 
larization direction has been shown to be 

of the exponential form and is very close to P°J^{cos9) 
for efc in Figs.7 and 8. Here, rik and are nearly par¬ 
allel from fj-oEo ^ SOfeeT. We have also calculated the 
distribution for the angle (p = cos“^(efc • n^) between rifc 
and efc, which indeed has a sharp peak at (^ = 0 in Fig.9. 
In addition, in Sec.HIC, we have shown that the relative 
vector — rkHi — '>'kH 2 between the two protons tends 
to be parallel to its conjugate force in Eq.(70). 

(5) In Sec.HID, we have calculated the single-molecule 
orientation time-correlation functions in Figs. 11 and 
12. Those for ek{t) and nk{t) decay exponentially as 
exp(—t/rh) with = 5.0 ps for t > 1 ps, which is gov¬ 
erned by the hydrogen bond dynamics. However, the 
perpendicular part of rikit) with respect to e.k(t)) under¬ 
goes librational motions on a timescale of 0.01 ps. We 
have also calculated those related to the relative vector 
^k(0 between the two protons to find similar results. 

We make some critical remarks, i) We should ex¬ 
amine how the dipole and the local field move to¬ 
gether in rotational big jumps in liquid wateil^TliSl, jj) 
We should consider the molecular polari zabil ity due to 
molecular stretching in strong local fielcPi^. hi) Ions 
strongly influence the local fields of the surrounding wa¬ 
ter molecule^^, resulting in a hydration shell around each 
ion, so this aspect should further be studied, iv) In fu¬ 
ture work, we will show that the dipoles Hk(t) are more 
aligned along the local field F)°‘^(t) in supercooled states 
than at T = 298 K. Thus, the slowing-down of the orien¬ 
tational dynamics in supercool ed wa ter stems from that 
of the hydrogen bond dynamic^^^^. 
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Appendix A: Local electric field and polarization 
density for water 

We explain how the local field and the polarization 
density are defined. In the TIP4P/2005 modeP^, each 
water molecule k forms a rigid isosceles triangle, where 
aoH = \rkm - ffcol = kfeH 2 - ^fcol = 0.957 A with the 
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HOH angle being 0 hoh = 104.5°. The charge position 
M is given by r^M = J’fcO with oom = 0.1546A. 

The rifc is the unit vector from r^o to the midpoint of 
the protons. The dipole moment is written as Eq.(l). 

The potential in Eq.(4) depends on the center of 
mass VkG = ^rko + + T-fcH 2 ), the relative vector 

between the two protons = '>'/cHi~'>'feH 2 j and the dipole 
moment For small shifts of the charge positions 
at hxed E^,, the potential change dU^n is rewritten ad^ 

du^ = -J2[FI ■ drkG + n • d^k + ■ dt^k]- (Al) 

k 

The conjugate electric force to VkG is given by = 
qH_{Ekkii + FfeH 2 - 2FfcM) and that to is given in 
Eq.(70). Conjugate to is the local electric held, 

A -®feH2) — by[EkM, (A2) 

where 6 m is a small coefficient given by 

6m = Soom/Somh ~ i/9 “ 0.208. (AS) 

Here, Omh = non cos(0hoh/2) — oom is the distance be¬ 
tween the point M and the midpoint of the protons, while 
1/9 is the protons-to-molecule mass ratio. 

The polarization density p{r) in Eq.(19) consists of the 
contributions from the constituting molecules as p(r) = 
J2kPkif')- For the TIP4P/2005 model, we havd^ 

Pkir) = ^[6(r;rfcHi,rfcH) - 5(r; rfcH2, rfcH)]gHlfc 

+S{r-rkH,rkM)Pk> (A4) 

where = (’’feHi + ’'fcH 2)/2 is the midpoint of the 
proton positions, and S{r-ri,r 2 ) = Jq dA S(r — Xri — 
(1 — A)r 2 ) is the symmetrized 6 functiorP^. We conhrm 
Eq.(19) and obtain the dipole and quadrupole moment 
as f drpf^ = /Xfc and J dr{r - rM)Pk = + 

4 Q}i^k^k ■ 

Appendix B: Polarization fluctnations in the 
presence of parallel metal walls 

We examine the polarization huctuations in highly po¬ 
lar huids without ions between metal walls. We use the 
continuum electrostatics, so the polarization density p 
and the electric held E are smooth variables in the bulk 
(outside the Stern layers). We assume that the cell length 
H much exceeds the Stern layer thickness d for simplicity. 
Then, the volume L^{H — 2d) of the bulk region is close 
to the cell volume V = L^H. In terms of % in Eq.(24) 
and C in Eq.(28), the electrostatic energy is written as 


{0 < x,y < L). The hrst term is the bulk contributiorP, 
the second term is due to the potential drop in the Stern 
layer^^, and the third term arises from the hxed poten¬ 
tial condition. In equilibrium, in Eq.(Bl) is mini¬ 
mized for p = Ph^z and E = E^Ez in the bulk region 
with ((To) = eFb/47r (see Eq.(21)). 

In this appendix, we superimpose deviations homoge¬ 
neous in the xy plane on the above equilibrium averages. 
That is, we pick up the Fourier components of the devi¬ 
ations with zero wave vector {kx,ky) = (0,0) in the xy 
plane. Then, Spz{z) = Pz — Ph depends only on z. As the 
electric held, we consider the huctuating Poisson electric 
held along the z axis, whose deviation is given by 

6Ez{z) = 47r(5(to — AnSpziz). (B2) 

Here, the deviation of the electric induction 6Dz{z) = 
SEz{z) + 47rSpz(z) is independent of z (in the absence 
of ions). The relation 6pz{z) = xdEz{z) needs not hold. 
The hxed potential condition yields the bulk averages, 

lEz = [ dz6Ez{z)/H =-6do/CH, (B3) 

Jh 

Spz = [ dz6pz{z)/H = {1 + iy,/eH)Sdo, (B4) 

Jh 

where 5p^ is the average of the deviation in the bulk, 

dPz = [ dr{pz{r) - Ph)/V. (B5) 

Jh 

We then calculate the change in in Eq.(Bl) in the 
bilinear (lowest) order in these deviations. The bulk av¬ 
erages in Eqs.(B3) and (B4) (oc Sd^) give 

6C/F) = ^(^hf + ^(^b)" + 

= B(1 + + l^/eH){5do)^l2x 

= V(1 + C/F6)(l + i^/eH)-\S^jy2x. (B6) 

In the hrst line, the bulk volume is equated with V. In the 
second and third lines, l^jePl ^ alPl ^ 1, so this factor 
can be omitted. Namely, we can simply set 5p^ = Sag and 
SEz = 0 in Eq.(B6). Next, the inhomogeneous deviations 
dpz — dpz: and 5Ez — SEz = —4Tr{Spz — 5p^) yield 

=L^ [ dz^iSpz-S^zY- (B7) 

Jh 2x 

The total change in is the sum SUm = SUm^ + SUm^ 
and the deviation 5pz{z) obeys the Gaussian distribu¬ 
tion cx exp(—(JGm/fcBT). We can then calculate the pair 
correlation function in the bulk along the z axis: 





\Er 

Stt 


dr_L^-L^A4>cro, (Bl) 


A(zi,Z2) 


dru_ 


dr2±Gzz{ri,r2)/kBTL'^, (B8) 


where J^dr is the integral in the bulk region {d < z < where G^^(ri,r 2 ) is the zz component of the 3D polar- 
H — d) and f dr± = f dxdy is the surface integral ization correlation function in Eq.(30) with zi and Z 2 in 
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the bulk region. From Eqs.(B6) and (B7), the ID corre¬ 
lation function K{zi,Z 2 ) satisfies the integral equation, 

-[X(zi,Z 2 ) -K] + -(1 + = S(zi - Z2). (B9) 

X X ^ 

Here, K = ^^dz\K{zi, Z 2 )/H is the average over z\ in¬ 
dependent of Z 2 (see Eq.(B12)). We solve Eq.(B9) as 


K{zi,Z2) = - 5 {zi - Z2) -f 

e ti -\- 


X 

eH’ 


(BIO) 


From Eq.(B12) below, the above relation is consistent 
with Eqs.(33) and (36) with = 1. In fact, the sum 
of the last two terms in Eq.(BlO) is ho/H with ho being 
given by Eq.(39). We have also obtained the coefficient 
x/e in front of 6 {zi—Z 2 )- This is because we have treated 
the dipolar interaction as a short-range one (along the 2 ; 
axis) in setting up Eq.(Bl). 

In the Stern layers, the microscopic polarization varies 
noticeably on the scale of d (see Fig.2). The integral of 
its 2 ; component in the layers is related to C bjim 


'stern 47rG 


tJo = - 


eH 


o'o, 


(BlI) 


which is equal to drp^ and is consistent with 

Eq.(B4). Since i^/eH «C 1, the integral of 5pz in the 
Stern layers is much smaller than that in the bulk region. 
Thus, we have drSpz = 5Mz, leading to 


/ dziK{zi,Z2) 

Jh 


ksTV 


X 

I+fw/H’ 


(B12) 


which confirms the consistency of Eqs.(23), (39), and 
(BIO). 


Appendix C: Dielectric formula for a polar fluid 
in a sphere 


In the continuum electrostatics, we consider the po¬ 
larization fluctuations in a polar fluid in a sphere with 
radius R much longer th an cr , which is embedded in an 
infinite dielectric mediunPE^. For the sake of generality, 
we assume that the dielectric constant of the sphere ex¬ 
terior e' can be different from that in the interior e. If 
the electric field tends to along the z axis far from 
the sphere, the average electric potential (j) is given by 


<(> 

E\yZ 2£‘^ “h c 


9{R — r) + 9{r — 



(Cl) 


where pin and Pex represent the fluctuation strength of 
the polarization in the sphere interior and exterior, re¬ 
spectively. Note that dcj) is continuous at r = i? in 
Eq.(C4). We further require the continuity of the elec¬ 
tric induction 5D = SE + Airdp along the surface normal 
(oc r) at r = R, where SE = —V5(f>. Then, the strength 
of the electric field deviation is determined as 


Ein = -47r(pi„ -f 2pex)/3. (C5) 


We use the electrostatic energy U^a in Eq.(Bl) neglect¬ 
ing the surface capacitance term. In the bilinear order of 
the deviations, its incremental change is written as 


SU^ = / dr 


-^\Sp\^ + Ese\^ 

2xe{r) Stt 


(C6) 


From V • SD = 0 we have J drSE ■ SD — 0 and 


SUm = / dr 


2Xe(r) 


6p-[6p-Xeir)SE]. (C7) 


Substitution of Eqs.(C3) and (C4) gives 


SUm 

V 


2e' + e 
2e' -f 1 


pI I + 1 
2x 3x' 


e'-l 


(C8) 


where v = AttR?/‘S. We recognize that pin and the 
combination q = pex +Pin(£' — l)/(2e' + 1) are indepen¬ 
dent Gaussian random variables obeying the distribution 
(X ejq){—SUni/k^T). In the linear response regime, we 
can set drSpz = vpin to obtain 

{{SMif) = vkBTx{2e' + l)/{2e' + e), (C9) 

which is well-known in the literatur (j^ l ^°b4 [ jf .^g 

reproduce Eq.(40) for v <^V. 

Using simple arguments, FrohlichP^ obtained the first 
term in Eq.(C8) for e' = e. To reproduce his result, let 
us set g = 0 at fixed pin. Then, Ei^ = —47rpin/(2e' +1) = 
Pexlx' from Eq.(C5) so that Sp—XeSE vanishes for r > i? 
and is equal to [(e + 2e')l{2e' + l)]pine-z for r < R. Thus, 
the outer region is in a local equilibrium stat e for q = 0 , 
which was assumed in the reaction field theorjPUH. Now, 
from Eq.(C7), the first term in Eq.(C8) readily follows. 

Appendix D: Polarization correlation function 


The average polarization is given by p = —Xe{'<’)XJ(j) in 
terms of the r-dependent susceptibility, 

Xe{r)=xd[R-'r)+x!d{T'- R), (C2) 


Here, we calculate the self part of the polarization 
correlation function in Eq.(30). In terms of Pfc(r) in 
Eq.(A4), it is expressed as 


where x' = ~ l)/47r. 

On the above equilibrium profiles, we superimpose 
small fluctuations of the polarization and the electric po¬ 
tential as in Appendix B. They are expressed as 

Sp = pioO{R -r)ez+ Po^9{r - R)R^V ^, (C3) 

^ = -E;oz[ 6 {R-r) + 6 {r-R)R^^l (C4) 


Gcp{ri,r2) = ^{Pka{ri)pip{r2)) 
k,i 

= G““(t’i.’'2)+ GP“(t'i,t-2), (D1) 

where we assume Ea = 0. In the second line, we divide 
Ga/j into the self part G^^^ with k = i and the pair part 

k^ 
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FIG. 15: (a) Functions ipir) and C(r) in Eqs.(D4) and (D5) 

vs. r, where 2A = ohu (HH distance) and h = ajyjg (distance 
between M and the proton midpoint), (b) Radius-dependent 
correlation factor Gk{R) in Eq.(D6) vs. R, which tends to 
the correlation factor for R > 2 K. It has a majdmum at 
r = 2A/3 due to the unique shape of the water molecules. 


For point-dipole models, we would have ^ 2 ) = 

— r2)/3. For the TIP4P/2005 modef^, the 
self part depends on the following molecular lengths, 

h = omh = 0.431 A, A = aHH/2 = 0.757 A, (D2) 


where h is the distance between the point M and the 
midpoint of the two protons, and A is half of the dis¬ 
tance between the two protons. The averages over the 
orientations of rik and give 




nfil r/ 3C\ XaX/S _ Sap 

2'Kh‘^ \ Ar J r‘^ 4r^ 


(D3) 


where r = r 2 — ri. We define = ip{r) and ^ = ({r) as 

'0(r) = (h — r)9{h — r) + {X — r)6{X — r) 

-(2A-r)6»(2A-r)/4, (D4) 

^(r) = -{h^ - r^)e{h - r) - (A^ - r‘^)9{X - r) 

+ {e^ — r‘^)9{e — r), (D5) 


where 9{u) is the step function being equal to 1 for it > 0 
and 0 for u < 0 and e = = 0.871 A. The self 

part is nonvanishing only for r < 2A = 1.51 A, where 
the upper bound is the H-H distance ohh- As shown in 
Fig.15(a), 'ip{r) > 0 for r < 2A/3 and ip{r) < 0 for 2A/3 < 
r < 2A leading to dr'il){r) = h?/2, while (^{r) > 0 for 
any r. Therefore, the integral of G^®^^(ri,ri -|- r) with 
respect to r is equal to n/Xg(5a/3/3, as ought to be the case. 

Following the previous author calculate the 
radius-dependent correlation factor defined by 


Gk{R)=[ dr'^Gaa{ri,ri+r)/nfxl. 

Jr<R 


(D6) 


In Fig.15(b), we plot Gk{R) vs the radius R using the 
self part in Eq.(D3) and a point-dipole approximation 
for the pair parlP^. We confirm the plateau behavior 
GK(i?) =gK = 2.14 for cr < i? < pi/3. 
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